CFA como modelo mixto

Ajustando un CFA para obtener resultados de un MLM

Author

https://dacarras.github.io/

Modified

November 14, 2024

Resumen

Simulamos datos multinivel con componentes de varianza de \(\sigma^2_{w} = .1\) y \(\sigma^2_{b} = .9\), y con media cero. Generamos datos para 500 clusters, con 5 observaciones para cada cluster, dando un total de 25000 observaciones. Estos datos, se pueden bajar del siguiente link.
Cada cluster lo simbolizamos con la letra j, y a cada observación con la letra i.
Los datos generados se encuentran en formato stacked o long, donde cada fila es una observación, indexada en un cluster. Estos datos los volteamos, para producir una representación de las mismas observaciones en formato wide o hacia el lado. En este caso, se emplean 5 columnas (i1-i5) para alojar los valores de cada observación por cada cluter j.
Empleando a la libreria lavaan, ilustramos como ajustar un modelo MLM sobre los datos en stacked, y ajustamos un modelo CFA sobre los datos en wide, recuperando los mismos resultados en cada caso, hasta a dos decimales.

Flujo de trabajo

1. Datos simulados
1. Preparar los datos
1. Ajustar Modelos
1. Tabla comparada de resultados

1. Datos Simulados

\[y_{ij} = \gamma_{00} + \mu_{0j} + \epsilon_{ij}\]

Donde,

\(y_{ij}\) = respuesta \(y\) para la observación \(_{i}\), en el cluster \(_{j}\)
\(\gamma_{00}\) = media entre todos los clusters
\(\mu_{0j}\) = desviaciones de las medias de los clusters, a la gran media \(\gamma_{00}\)
\(\epsilon_{ij}\) = desviación de las observaciones con respecto a la media de su propio cluster
\(\mu_{0j} \sim \mathcal{N}(0,\,\sigma^{2}_{b})\) = donde \(\mu_{0j}\) es variable de media cero, y varianza \(\sigma^{2}_{b}\)
\(\sigma^{2}_{b}\) = .9
\(\epsilon_{0j} \sim \mathcal{N}(0,\,\sigma^{2}_{w})\) = donde \(\epsilon_{ij}\) es variable de media cero, y varianza \(\sigma^{2}_{w}\)
\(\sigma^{2}_{w}\) = .1
\(_{j}\) = 500 clusters
\(_{i}\) = 5 observaciones para cada cluster

Datos simulados:

data_long.rds

Nota

Los datos simulados fueron generados empleando al software Mplus (ver Anexo).

2. Preparar datos

Abrimos los datos simulados, y generamos los datos en formato wide

# ubicación de los datos en url
url_file_long <- url('https://github.com/dacarras/edu4046_cfa/raw/main/data_long.rds')

# datos en formato stacked
data_long <- readRDS(url_file_long)

# datos en formato wide
data_wide <- data_long %>%
             mutate(id_ij = paste0('i',rep(seq(1:5),500))) %>%
             tidyr::pivot_wider(
              id_cols = -id_i,
              names_from = 'id_ij', 
              values_from = 'y') %>%
             dplyr::glimpse()

Rows: 500
Columns: 6
$ id_j <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19…
$ i1   <dbl> -0.649609, 1.810764, -1.365712, -0.555940, -1.024259, -0.441473, …
$ i2   <dbl> -0.513657, 0.345480, -1.018105, -0.355253, -1.083364, -0.886687, …
$ i3   <dbl> -0.697763, 0.818963, -1.361241, -0.655227, -1.675723, -0.422175, …
$ i4   <dbl> 0.009125, 0.982486, -1.054176, 0.053315, -1.880564, 0.001164, 2.1…
$ i5   <dbl> 0.344252, 0.929602, -0.563126, -0.030324, -1.281889, -1.597010, 1…

Cómo se ven los datos generados, y los datos en formato wide, para los 4 primeros clusters.

# datos en formato stacked
data_long %>%
dplyr::slice(1:20) %>%
knitr::kable(., digits = 4)

id_i	id_j	id_ij	y
1	1	1	-0.6496
2	1	2	-0.5137
3	1	3	-0.6978
4	1	4	0.0091
5	1	5	0.3443
6	2	1	1.8108
7	2	2	0.3455
8	2	3	0.8190
9	2	4	0.9825
10	2	5	0.9296
11	3	1	-1.3657
12	3	2	-1.0181
13	3	3	-1.3612
14	3	4	-1.0542
15	3	5	-0.5631
16	4	1	-0.5559
17	4	2	-0.3553
18	4	3	-0.6552
19	4	4	0.0533
20	4	5	-0.0303

# datos en formato wide
data_wide %>%
dplyr::slice(1:4) %>%
knitr::kable(., digits = 4)

id_j	i1	i2	i3	i4	i5
1	-0.6496	-0.5137	-0.6978	0.0091	0.3443
2	1.8108	0.3455	0.8190	0.9825	0.9296
3	-1.3657	-1.0181	-1.3612	-1.0542	-0.5631
4	-0.5559	-0.3553	-0.6552	0.0533	-0.0303

3. Ajustar modelos

3.1 Modelo MLM tradicional

Empleamos la función lme4::lmer() para ajustar un modelo multinivel de tipo nulo, sobre los datos data_long.rds.

# ----------------------------------------------- 
# method 1: lem4::lmer()
# -----------------------------------------------

fit_null <- lme4::lmer(y ~ 1 + (1|id_j), data = data_long, REML = FALSE)
summary(fit_null)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: y ~ 1 + (1 | id_j)
   Data: data_long

     AIC      BIC   logLik deviance df.resid 
  3315.2   3332.6  -1654.6   3309.2     2497 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0237 -0.6097 -0.0063  0.6006  3.7461 

Random effects:
 Groups   Name        Variance Std.Dev.
 id_j     (Intercept) 0.9564   0.9780  
 Residual             0.1013   0.3183  
Number of obs: 2500, groups:  id_j, 500

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.006864   0.044196   0.155

3.2 Modelo MLM en lavaan

Empleamos la función lavaan::sem() para ajustar un modelo multinivel de tipo nulo, sobre los datos data_long.rds.

# ----------------------------------------------- 
# method 2: lavaan::sem()
# -----------------------------------------------

lavaan_mlm <- '
level: 1
# variance w
y~~w*y

level: 2
# variance b
y~~b*y

y~mu*1

# model terms
y00   := mu
var_w := w
var_b := b
icc   := b/(b+w)

'

fit_mlm <- lavaan::sem(model = lavaan_mlm, 
           data = data_long, 
           cluster = "id_j"
           )

summary(fit_mlm)

Length  Class   Mode 
     1 lavaan     S4

lavaan::parameterEstimates(fit_mlm) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
dplyr::select(label, est, se) %>%
knitr::kable(., digits = 3)

label	est	se
y00	0.007	0.044
var_w	0.101	0.003
var_b	0.956	0.062
icc	0.904	0.006

3.2 Modelo MLM por medio de un modelo CFA

Empleamos la función lavaan::sem() para ajustar un modelo multinivel de tipo nulo a través de un modelo CFA, sobre los datos data_wide.

# ----------------------------------------------- 
# method 3: lavaan::sem()
# -----------------------------------------------

lavaan_cfa <- '
# random term
eta =~ 1*i1 
eta =~ 1*i2 
eta =~ 1*i3 
eta =~ 1*i4 
eta =~ 1*i5 

# latent mean
eta ~ mu*1

# person intercepts
i1  ~ d1 *1
i2  ~ d2 *1
i3  ~ d3 *1
i4  ~ d4 *1
i5  ~ d5 *1

# variance
eta~~b*eta

# residual variance
i1 ~~ w*i1
i2 ~~ w*i2
i3 ~~ w*i3
i4 ~~ w*i4
i5 ~~ w*i5

# model constraints

cw := (d1 +d2 +d3 +d4 +d5)/5

cw == 0

# model terms
y00   := mu
var_w := w
var_b := b
icc   := b/(b+w)


'

fit_cfa <- lavaan::sem(model = lavaan_cfa, 
           data = data_wide
           )

summary(fit_cfa)

Length  Class   Mode 
     1 lavaan     S4

lavaan::parameterEstimates(fit_cfa) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
dplyr::select(label, est, se) %>%
knitr::kable(., digits = 3)

label	est	se
y00	0.007	0.044
var_w	0.101	0.003
var_b	0.956	0.062
icc	0.905	0.006

4. Tabla comparada de resultados

4.1 Extraer los estimados de cada modelo ajustado

# ----------------------------------------------- 
# retrive estimates
# -----------------------------------------------

# function to adjust to two decimal places
decimal <- function (x, k) {
    format(round(x, k), nsmall = k)
}

# lmer4::lmer()
table_est_1 <- broom.mixed::tidy(fit_null) %>%
mutate(label = c('y00','var_b', 'var_w')) %>%
mutate(est = case_when(
effect == 'fixed'    ~ estimate,
effect == 'ran_pars' ~ estimate*estimate
)) %>%
mutate(est_1_se = paste0(decimal(est, 2), ' (', decimal(std.error,2), ')')) %>%
dplyr::select(label, est_1_se)

# lavaan mlm
table_est_2 <- lavaan::parameterEstimates(fit_mlm) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
mutate(est_2_se = paste0(decimal(est, 2), ' (', decimal(se,2), ')')) %>%
dplyr::select(label, est_2_se) 

# lavaan cfa
table_est_3 <- lavaan::parameterEstimates(fit_cfa) %>%
dplyr::filter(op == ':=') %>%
dplyr::filter(label %in% c('y00','var_w','var_b','icc')) %>%
mutate(est_3_se = paste0(decimal(est, 2), ' (', decimal(se,2), ')')) %>%
dplyr::select(label, est_3_se)

4.2 Tabla comparada

# ----------------------------------------------- 
# display table comparison
# -----------------------------------------------

dplyr::full_join(table_est_1, table_est_2, by = 'label') %>%
dplyr::full_join(., table_est_3, by = 'label') %>%
rename(
`lme4::lmer()` = 2,
`lavaan::sem()-MLM` = 3,
`lavaan::sem()-CFA` = 4
) %>%
knitr::kable()

label	lme4::lmer()	lavaan::sem()-MLM	lavaan::sem()-CFA
y00	0.01 (0.04)	0.01 (0.04)	0.01 (0.04)
var_b	0.96 ( NA)	0.96 (0.06)	0.96 (0.06)
var_w	0.10 ( NA)	0.10 (0.00)	0.10 (0.00)
icc		0.90 (0.01)	0.90 (0.01)

Anexos

Código Mplus para generar los datos simulados

TITLE:
null mlm;

MONTECARLO:
 
NAMES         = y;       
NOBSERVATIONS = 2500;    
NCSIZES       = 1;       
CSIZES        = 500 (5); 
NREPS         = 1;       
SEED          = 4046;    
SAVE          = icc_90.dat;

MODEL POPULATION:
 
%WITHIN%                 
!residual variance       
y*0.1;

%BETWEEN%               
!latent mean            
[y@0];                  
y*0.9;

ANALYSIS:

TYPE = TWOLEVEL;
ESTIMATOR = ML;

MODEL:
 
%WITHIN%                 
!residual variance       
y*0.1;

%BETWEEN%               
!latent mean            
[y@0];                  
y*0.9;

OUTPUT:
TECH9;